function [A, b] = process_boundary_2d(A, b, mesh, pde)

ref_gauss_points_1d = [-0.9681602395076260898355762, -0.8360311073266357942994298, -0.6133714327005903973087020, -0.3242534234038089290385380, 0, 0.3242534234038089290385380, 0.6133714327005903973087020, 0.8360311073266357942994298, 0.9681602395076260898355762];
ref_gauss_weights_1d = [0.0812743883615744119718922, 0.1806481606948574040584720, 0.2606106964029354623187429, 0.3123470770400028400686304, 0.3302393550012597631645251, 0.3123470770400028400686304, 0.2606106964029354623187429, 0.1806481606948574040584720, 0.0812743883615744119718922];
psi = reference_basis_2d(mesh.basis_type);

% grad_u = [pde.dudx; pde.dudy];
% p = grad_u'*vector_normal;
% 1, vector_normal = [0; -1];
% 2, vector_normal = [1; 0];
% 3, vector_normal = [0; 1];
% 4, vector_normal = [-1; 0];
%% Neuman
vector_b = zeros(mesh.Nb,1);
for k = 1:size(pde.bdary, 2)
    if pde.bdary(1,k) == 1
        bde = mesh2boundary_edges(mesh, k);
        switch k
            case 1
                for i = 1:size(bde,2)
                    nk = bde(1,i);
                    E = mesh.node(mesh.elem(nk,:),:);
                    V = E(1:3,1:2);
                    J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
                    x_loc = @(x,y) ((x-V(1,1))*(V(3,2)-V(1,2)) - (y-V(1,2))*(V(3,1)-V(1,1)))/J_det;
                    y_loc = @(x,y) ((V(2,1)-V(1,1))*(y-V(1,2)) - (V(2,2)-V(1,2))*(x-V(1,1)))/J_det;
                    p1 = mesh.node(bde(2,i),1);
                    p2 = mesh.node(bde(3,i),1);
                    gauss_points_1d = 0.5*(p2-p1)*ref_gauss_points_1d + 0.5*(p1+p2);
                    gauss_weights_1d = 0.5*(p2-p1)*ref_gauss_weights_1d;
                    for beta = 1:size(mesh.elem,2)
                        r = 0;
                        for j = 1:9
                            r = r + gauss_weights_1d(j) * pde.cor(gauss_points_1d(j), mesh.bottom) * (-1) * pde.dudy(gauss_points_1d(j), mesh.bottom) * psi{beta}(x_loc(gauss_points_1d(j), mesh.bottom), y_loc(gauss_points_1d(j), mesh.bottom));
                        end
                        vector_b(mesh.elem(nk,beta),1) = vector_b(mesh.elem(nk,beta),1) + r;
                    end
                end
            case 2
            case 3
            case 4
        end
    end
end
b = b + vector_b;

%% Dirichlet
for k = 1:size(pde.bdary, 2)
    if pde.bdary(1,k) == 0
        bdn = mesh2boundary_nodes(mesh, k);
        A(bdn,:) = 0;
        A(sub2ind(size(A), bdn, bdn)) = 1;
        b(bdn) = pde.u(mesh.node(bdn,1),mesh.node(bdn,2));
    end
end

end

%% boundary_nodes
function bdn = mesh2boundary_nodes(mesh, index)
bd_bottom = 1:mesh.Ny+1:mesh.Np;
bd_right = mesh.Np-mesh.Ny:mesh.Np;
bd_top = mesh.Np:-(mesh.Ny+1):mesh.Ny+1;
bd_left = mesh.Ny+1:-1:1;
switch mesh.basis_type
    case {"P2", "Q2"}
        bd_bottom = [bd_bottom, mesh.Np+mesh.Ny+1:3*mesh.Ny+1:mesh.Nb];
        bd_right = [bd_right, mesh.Nb-mesh.Ny+1:mesh.Nb];
        bd_top = [bd_top, mesh.Nb-mesh.Ny:-(3*mesh.Ny+1):mesh.Np+1];
        bd_left = [bd_left, mesh.Np+mesh.Ny:-1:mesh.Np+1];
end

switch index
    case 1
        bdn = bd_bottom;
    case 2
        bdn = bd_right;
    case 3
        bdn = bd_top;
    case 4
        bdn = bd_left;
    otherwise
        error('Invalid edge number');
end

end

%% boundary_edges
function bde = mesh2boundary_edges(mesh, index)
bd_bottom = [1:2*mesh.Ny:mesh.Ne;
    1:mesh.Ny+1:mesh.Np-mesh.Ny-1;
    mesh.Ny+2:mesh.Ny+1:mesh.Np];
bd_right = mesh.Np-mesh.Ny:mesh.Np;
bd_top = mesh.Np:-(mesh.Ny+1):mesh.Ny+1;
bd_left = mesh.Ny+1:-1:1;

switch index
    case 1
        bde = bd_bottom;
    case 2
        bde = bd_right;
    case 3
        bde = bd_top;
    case 4
        bde = bd_left;
    otherwise
        error('Invalid edge number');
end

end